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Abstract —We develop a novel sampling theorem for functions 
defined on the three-dimensional rotation group SO(3) hy con¬ 
necting the rotation group to the three-torus through a periodic 
extension. Our sampling theorem requires 4L® samples to capture 
all of the information content of a signal band-limited at L, 
reducing the number of required samples by a factor of two 
compared to other equiangular sampling theorems. We present 
fast algorithms to compute the associated Fourier transform on 
the rotation group, the so-called Wigner transform, which scale as 
0(L'^), compared to the naive scaling of 0{L^). For the common 
case of a low directional band-limit N, complexity is reduced to 
0{NL^). Our fast algorithms will be of direct use in speeding up 
the computation of directional wavelet transforms on the sphere. 
We make our S03 code implementing these algorithms publicly 
available. 

Index Terms —Harmonic analysis, sampling, spheres, rotation 
group, Wigner transform. 

I. Introduction 

HANNON established the theoretical foundations of sam¬ 
pling theory in Euclidean space over half a century ago, 
proving that the information content of a band-limited contin¬ 
uous signal could be captured completely in a hnite number 
of samples [1]. The fast Fourier transform (EFT) [2] is one of 
the most important algorithmic developments of our era and 
has been instrumental in rendering the frequency content of 
signals accessible in practice. The combination of theoretical 
foundations and fast algorithms has led to the extensive use 
of Fourier methods to analyse data in myriad applications. 

In many applications, however, data are acquired on non- 
Euclidean manifolds where Euclidean sampling theory is not 
applicable. Spherical manifolds are one of the most prevalent 
non-Euclidean domains. When observing over directions, data 
are acquired on the two-dimensional sphere. If distance infor¬ 
mation is also accessible, then data are acquired on the three- 
dimensional ball. For example, in cosmology observations are 
made on the celestial sphere, with the cosmic microwave 
background (CMB), observed on the sphere {e.g. [3]), while 
surveys of the distribution of galaxies are made on the ball {e.g. 
[4]). Other examples of data dehned on spherical manifolds 
are found in many helds, including planetary science {e.g. [5], 
[6]), geophysics {e.g. [7]) and biomedical imaging {e.g. [8]). 
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Sampling theory and corresponding fast harmonic trans¬ 
forms on the sphere remain active areas of research. Until 
recently, the canonical equiangular sampling theorem on the 
sphere was that of Driscoll & Healy [9], requiring ^ 4L^ 
samples on the sphere to capture the information content of 
a signal band-limited at L. A novel sampling theorem on the 
sphere was developed recently by two of the authors of the 
present article (McEwen & Wiaux), reducing the number of 
samples required to capture a band-limited signal to ~ 2L^ 
[10], building on the developments of [11], [12] (see [10] 
for a review of sampling theory on the sphere). No existing 
sampling theorem on the sphere reaches the optimal number 
of samples given by the harmonic dimensionality of a band- 
limited signal of L^. A new sampling scheme that achieves 
the optimal number of samples was developed recently by 
[13]. Whilst this scheme does not lead to a sampling theorem 
with theoretically exact spherical harmonic transforms, good 
numerical accuracy is achieved in practice and fast algorithms 
were developed. A new sampling theory on the ball was 
developed by [14] recently, augmenting the sampling theorem 
on the sphere of [10] with Gaussian quadrature on the radial 
line to recover an exact harmonic transform suitable for the 
analysis of large data-sets dehned on the ball. 

The analysis of data dehned on the sphere or ball often 
leads to data dehned on the rotation group SO(3), the space 
of three-dimensional rotations. Directional wavelet transforms 
on the sphere {e.g. [15]-[24]) probe signal content not only in 
scale and position on the sphere, but also in orientation. The 
resulting wavelet coefficients thus live on the rotation group. 
Consequently, in all applications where directional spherical 
wavelet transforms are performed, sampling on the rotation 
group is important. Moreover, the wavelet transform can be 
computed via a Fourier transform on the rotation group [21], 
[22]. Data dehned natively on the rotation group also arise 
in many applications, for example searching databases of 
objects over arbitrary rotations [25]. Sampling theorems on 
the rotation group with fast harmonic transforms are thus of 
both important theoretical interest and practical use. 

The canonical equiangular sampling theorem on the rotation 
group SO(3) is that of Kostelec et al. [26], which relies on 
the Driscoll & Healy sampling theorem on the sphere [9], and 
requires ^ 8L^ samples to capture a signal on the rotation 
group that is band-limited at L. In this article we develop a 
novel sampling theorem on the rotation group, reducing the 
number of samples required to capture a band-limited signal 
to ~ 4L^. Furthermore, we present fast algorithms to compute 
the associated harmonic transform on the rotation group (often 
called the Wigner transform). No existing sampling theorem 
on the rotation group reaches the optimal number of samples 
given by the ^ 4L^/3 degrees of freedom in harmonic space, 
although our approach comes closest to this bound. 
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The remainder of this article is structured as follows. After 
concisely reviewing harmonic analysis on the rotation group 
we present our novel sampling theorem, before defining ex¬ 
plicitly an exact quadrature rule on the rotation group. We then 
discuss and evaluate fast algorithms for computing the Wigner 
transform associated with our sampling theorem, presenting a 
new public code, followed by concluding remarks. 

II. Harmonic Analysis on the Rotation Group 

We consider the space of square integrable functions 
on the rotation group L^(SO(3)), with inner product 
if, g) = Iso{3) Mp) fip) 9*ip) for f,g G L2(S0(3)), 
where dg{p) = sin /3 da d/3 d 7 is the usual invariant mea¬ 
sure on SO(3), which is parameterised by the Euler angles 
p = (a,/ 3 , 7 ) G SO(3), with a G [0, 27r), /3 G [0,7r] and 
7 G [0, 27r). We adopt the zyz Euler convention corresponding 
to the rotation of a physical body in d. fixed coordinate system 
about the z, y and z axes by 7 , /3 and a, respectively. 

The Wigner 13-functions with natural f G N 

and integer m,n G Z, |m|, |n| < i, are the matrix 

elements of the irreducible unitary representation of the 
rotation group SO(3). Consequently, the also form 

an orthogonal basis in L^(SO(3)).' The orthogonality and 
completeness relations for the Wigner Z3-functions read, 
respectively, 

^£i' ^mm' ^nn' /{2i + 1) 

and ET=oT,L=-iEL-iDLn{c^,P,l)DtJa\fi',y) = 
i5(a — Q!')i5(cos/ 3 — cos/3') d (7 — 7 '), where dij is the Kro- 
necker delta symbol and S{x) is the Dirac delta function. The 
Wigner functions may be decomposed as [27] 

where the real polar d-functions can be computed by recursion 
{e.g. [28], [29]). The Wigner d-functions have the following 
Eourier series decomposition [30]: 

i 

dL„(/3) = i”-"* E ( 2 ) 

m' ——£ 

where = d^„( 7 r/ 2 ). This expression follows from a 

factoring of rotations as highlighted by Risbo [28]. We also 
note that the Wigner functions are related to the spin spherical 
harmonics by [31] 

/o/ -I- 1 

sYe^{9,y^) = {-ly ^. (3) 

A square integrable function defined on the rotation group 
may thus be represented by its Eourier expansion 

f(p) = E E E fLDtniP), (4) 

i—O m——in——£ 

where the Eourier coefficients are given by 

fL = (/, Dtj = f dgip) fip) Dl^ip) . (5) 

dS0(3) 

The Eourier transform on the rotation group SO(3) is often 
called the Wigner transform, while the Eourier coefficients 
are often called Wigner coefficients. 

'We adopt the conjugate D-functions as basis elements since this conven¬ 
tion simplifies connections to wavelet transforms on the sphere. 


III. Sampling Theorem 


We develop a novel sampling theorem on the rotation group 
by making a connecting to the three-torus,^ extending similar 
approaches on the sphere [ 10 ]-[ 12 ] (following closely the 
approach of our sampling theorem on the sphere [10]).^ Band- 
limited signals / G L^(SO(3)) are considered, with = 0 
for all f > L, TO > M and n > N (with M, N < L). Our 
sampling theorem is encapsulated in an exact computation of 
the Wigner transform of / from a finite set of samples. 

We adopt an equiangular sampling of the rotation group 
with sample positions given by 

where a G {0,1,..., 2M — 2} , ( 6 ) 


— 


2M- 1’ 


Pb = ’ where 6 g {0,1,...,L-1}, (7) 

and 

7g = 2 ^- 1 ’ where 5 G {0,1,... ,2A-2} , ( 8 ) 

(see [33] for alternative samplings of SO(3)). To make the 
connection with the three-torus, the domain is extended 
to [0,27r) by simply extending the domain of the b index 
to include {L, L -f 1,..., 2L — 1}. The number of required 
samples is thus [[L - 1)(2M - 1) + l] {2N - 1) ~ ALMN, 
or ~ samples when L = M = N. 

By noting the Wigner decomposition of Eqn. (1) and the 
Eourier series representation of Eqn. (2), the forward Wigner 
transform of Eqn. (5) may be written 


L-l 


y = i" 

J mn 


E a^ c 


^m'm ^m'n ^mnm' 5 


'=-(L-l) 


where 




/ dg{p) fip) Q-i(ma+m'p+nj) 
/sO(3) 

[ d/3sin/3G^„(/3)e-‘'”'^ 


(9) 

( 10 ) 

( 11 ) 


and 


p27T p2TT 

Gmniy) = da djfip) 
Jo Jo 


0 

M-1 


je 

N-l 


— i(ma+n 7 ) 


= (271)^ E E 

a=-(M-l)g=-(N-l) 

X e-‘l”'““+"^»V[(2M- 1)(2A- 1)] 


( 12 ) 


(13) 


Since Wigner coefficients are not defined for |m|,|n| > j 
we set them to zero to enforce the constraints |m|,|n| < 


^Although SO (3) is not homeomorphic to the three-torus T^, by excluding 
a subset of measure zero from SO(3) the resulting space can be embedded 
in T^, which is sufficient for sampling purposes. 

^Gauss-Legendre quadrature may also be used to define an efficient 
sampling theorem on the rotation group by a similar extension of the 
approach outlined in [10] and is developed in [32], where corresponding fast 
Wigner transforms are constructed by a separation of variables. Although the 
asymptotic number of samples is 4L^ in both cases, the approach described 
in this article requires fewer samples than Gauss-Legendre quadrature, which 
for small band-limits can be significant. 
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when the order of summations and integrals are interchanged. 
The final expression of Eqn. (13) follows by appealing to the 
orthogonality of the complex exponentials. 

To recover a sampling theorem on the rotation group we de¬ 
velop an exact quadrature for Eqn. (11) by extending GmniP) 
to the domain [0, 2tt), through the construction 


GmniP) 


GmniP) , /3e[0,7r] 

(_!)-+-/3e(7r,2^) 


(cf. [10], [12]), which may be represented by its Eourier 
decomposition 

L-l 

GmnW) = {2TTf ^ Gmnm' . (15) 

Substituting Eqn. (15) into Eqn. (11), one recovers the exact 
quadrature 


L-l 

Gmnm' ~ ^ ^ wi^Tfl ^ ) ; ( 16 ) 

m’' — — {L — 1) 


where the weights are given by w{m') = fj d^sin/3 e'™ 
which can be evaluated analytically [10]. 

The inverse Wigner transform may be computed in an 
analogous manner, again noting Eqn. (1) and Eqn. (2). The 
Eourier coefficients of the extension of / to the three-torus 
are first computed by 


L-l 




^ AC AC fC 
0—9 ^m'm ‘~^m'n J mn ’ 


^=0 




(17) 


followed by the three-dimensional Eourier transform 

M-l AT-l L-l 

f{aa,Pb,lg) = E E E Fmnm' 

m——(M —1) n——(N—l) m' —— (L —1) 

^ ^iimaa+m'fy+njg) 


The samples of / computed over /3 G (tt, 27r) are discarded. 


IV. Exact Quadrature 

Our sampling theorem on the rotation group can be used to 
define an explicit quadrature rule for the integration of a band- 
limited function. Approximately a quarter of the number of 
samples of the sampling theorem are required for the quadra¬ 
ture rule since integrating a band-limited signal corresponds to 
computing /qq only (aliasing in higher order coefficients can 
be tolerated). The exact quadrature reads: 

1= [ dg{p)f{p) (19) 

7sO(3) 

M-l L-l N-1 

= E E E fi^a^Pb,7g)Q{M , (20) 

a—0 0 g—0 


where = 2'Ka/M and = 2TTh/N (noting the reduced 
domain of a and h). The number of samples required is 
[{L-l)M + l]N r-. LMN, or ~ samples when L = 
M = N. The quadrature weights are defined by 


qiPb) = 


(27r 


|2 r 


MN I 


KPb) + (1 — <5b,L-l) u(/32L-2-h) 


( 21 ) 


where [10] 

= E (22) 

V. East Algorithms 

Fast algorithms to compute the Wigner transforms associ¬ 
ated with our sampling theorem can be implemented through 
a separation of variables, as described in Sec. Ill, using 
FFTs throughout to compute Fourier transforms efficiently. 
Furthermore, Eqn. (16) can be computed more efficiently in 
Fourier space, allowing aliasing in terms |m'| > L — 1 so 
that only L samples are needed in (3 (as described in [10]). 
Rather than implement this algorithm from scratch, however, 
we instead make the connection to spin spherical harmonic 
transforms to leverage our existing SSHT"^ code [10]. 

The forward Wigner transform may be expressed as 

/n(a,/3) = ^^ E 7s) (23) 

g=-{N-l) 

where 

I A 77 - /*27r 

V ani, ■''’""Vo 

x/„(a,/3)_„y;^(/3,a). (24) 

Eqn. (24) is simply a spin spherical harmonic transform 
with spin number —n, which may be computed for each n 
with asymptotic complexity 0{L^) by SSHT. The Fourier 
transform of Eqn. (23) can be computed by an EFT in 
0{L^N \ 0 g 2 N). The forward Wigner transform can thus be 
computed with overall complexity 0{NL^). 

The inverse Wigner transform may be expressed as 

L-l min(^,M-l) 
i=0 m=-min(r,M-l) 

(25) 

where 

N-l 

f{a,l3,jg)= fn{a,/3)F^^^ . (26) 

n=-(A-l) 

Eqn. (26) is simply an inverse spin spherical harmonic trans¬ 
form with spin number —n, which may be computed for 
each n with asymptotic complexity 0{L^) by SSHT. The 
Fourier transform of Eqn. (26) can be computed by an EFT 
in 0{Lp'N\og2 N). The inverse Wigner transform can thus be 
computed with overall complexity 0{NL^). 

For real functions on the rotation group we exploit the 
conjugate symmetry relation to 

reduce memory requirements and computational time by a 
factor of two. 

The fast Wigner transforms implementing our novel sam¬ 
pling theorem on the rotation group are implemented in our 
new S03^ code, which uses SSHT to compute spin spherical 
harmonics transforms and FFTW® to compute Fourier trans- 

"^http://www.spinsht.org 
^ http ://w w w. sothree.org 
®http://www.fftw.org 
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(a) N = L 



(b) AT = 4 

Fig. 1. Numerical accuracy of the algorithms implementing our sampling 
theorem for real and complex signals for N = L (panel a) and N = 4 
(panel b). Accuracy is found empirically to scale approximately linearly with 
increasing band-limit. 



(a) N = L 



(b) N = 4 

Fig. 2. Computation time of the algorithms implementing our sampling 
theorem for real and complex signals for N = L (panel a) and N = 4 (panel 
b). As predicted, computation time scales as O(NL^) and real transforms 
are approximately twice as fast as complex transforms (since additional 
symmetries are exploited). 


forms. The core code of S03 is written in C, while Matlab 
interfaces are also exposed. 

VI. Evaluation 

In order to assess the numerical accuracy and computation 
time of our algorithms we perform the following numerical 
experiment. We generate band-limited test signals defined by 
uniformly random Wigner coefficients with real and imaginary 
parts distributed in the interval [—1,1]. An inverse transform 
is performed to synthesise the test signal on the rotation group 
from its Wigner coefficients, followed by a forward transform 
to recompute Wigner coefficients. Numerical accuracy is mea¬ 
sured by the maximum absolute error between the original 
Wigner coefficients and the recomputed values. Computation 
time is measured by the average of the times taken to perform 
the inverse and forward transform. 

All numerical experiments are performed on a single core 
of a 2.67 GHz Intel(R) Xeon(R) CPU X5650 machine with 
100 GB of RAM and are averaged over ten random test 
signals. We consider two modes of operation: (i) N = L 
and (ii) a small constant N, in this example = 4. The 
latter mode of operation is a common use-case, for example 
when computing directional wavelet transforms on the sphere 
for wavelets with low azimuthal band-limit, and allows one 
to consider very large harmonic band-limits L. Note that 
the results presented here are obtained using version l.lbl 
of SSHT, which contains many computational optimisations 
compared to previous versions. The Risbo recursion [28] for 
computing Wigner c?-functions is adopted. 


The maximum absolute error is plotted against the band- 
limit in Fig. 1. High numerical accuracy is achieved, with 
errors on the order of machine precision and found empirically 
to increase approximately linearly with band-limit. 

The computation time for both real and complex signals 
on the rotation group is plotted against the band-limit in 
Fig. 2. Computation time evolves as 0{NL^) as predicted. 
Computation time for real signals is approximately twice as 
fast as for complex signals, also as predicted. 

VH. Conclusions 

We have presented a new sampling theorem on the rotation 
group SO(3) and fast algorithms to compute the associ¬ 
ated Wigner transform. Our sampling theorem requires 4L^ 
samples to capture the full information content of a signal 
band-limited at L, reducing the number of required samples 
by a factor of two compared to other equiangular sampling 
theorems on the rotation group [26]. Our fast algorithms are 
theoretically exact and achieve accuracy close to machine 
precision, while scaling as 0{L^), compared to the naive 
scaling of 0{L^). For the common case of a low directional 
band-limit N the scaling is reduced to 0{NL^). In a separate 
article [23] we apply these fast algorithms to improve the 
efficiency of a directional wavelet transform on the sphere 
[21], [22] in order to support the analysis of large data-sets. 
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